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The  Operator  Compact  Implicit  Method  for  Parabolic  Equations 


The  results  presented  in  this  report  were  obtained  as  part  of  a concerted 
effort  to  develop  reliable  and  efficient  numerical  techniques  for  solving 
major  fluid  dynamics  problems.  Here  the  methods  are  designed  to  form  the 
basis  for  a practicable  computer  code  to  solve  viscous  fluid  flow  problems. 
Efficient  fourth  order  finite  difference  approximations  are  developed.  Their 
associated  stability  and  accuracy  characteristics  are  analyzed  and  studied. 

These  techniques  presently  form  the  basis  for  a three-dimensional  boundary 
layer  computer  code  being  developed  at  NSWC/WOL. 
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T . INTRODUCTION 

The  current  engineering  requirements  for  providing  computational  fluid 
dynamics  codes  for  realistic  viscous  flow  problems  have  provided  the  impetus 
for  the  development  and  implementation  of  higher  order  finite  difference 
techniques  [61,  [ 1],  [22].  It  has  been  repeatedly  demonstrated  on  model 
problems,  that  even  the  simplest  types  of  higher  order  methods  should  provide 
tremendous  practical  advantages  in  terms  of  diminishing  the  required  number  of 
points  (storage)  and  also  the  overall  computing  time  for  a desired  resolution. 

The  present  effort  was  undertaken  to  confront  the  full  range  of  associated 
computational  problems  that  would  be  involved  in  practical  viscous  flow  field 
calculations.  Our  goal  was  to  try  to  develop  a cohesive  set  of  higher  order 
approximation  tools  which  would  help  to  indicate  what  methods  ultimately  might 
be  best  employed  to  form  the  basis  of  a major  new  code. 

It  appeared  to  several  people  almost  simultaneously  (sparked  by  a suggestion 
of  H.  0.  Kreiss  [14])  that  from  among  the  various  techniques  available  a fruitful 
class  of  methods  might  emerge  from  the  so-called  compact  implicit  techniques 
T 31,  [61.  Although  there  appear  to  be  a variety  of  forms  and  implementations 
the  approaches  do  share  some  broad  characteristics.  The  higher  order  is  usually 
sought  for  the  spatial  part  of  the  differential  operator.  The  method  developed 
is  generallv  required  to; 

1.  reduce  to  tridiagonal  form  for  fourth  order  accuracy 

2.  allow  for  nonuniform  spatial  grids  (usually  at  the  expense  of  one 
order  of  accuracy) 

3.  allow  for  flexibility  in  choosing  the  time  step. 

In  the  various  methods  developed  so  far  all  these  conditions  have  been  met  for 
simple  model  problems.  However,  further  important  concerns  still  remain. 

As  pointed  out  by  [ 3 ],  [ 4 ] and  [ 2 ] the  usual  compact  implicit  techniques, 
because  of  their  implicit  complexity,  are  not  generally  applicable  in  a direct 
manner  to  problems  with  varying  order  derivative  terms  unless  a vector  unknown 
of  the  derivative  values  is  considered.  Indeed,  adopting  the  factorization 
technique  suggested  in  [ 3 ] for  a wave  equation  problem  to  a model  parabolic 
problem  resulted  in  numerical  instabilities  (see  Section  III. 3 below).  To 
circumvent  such  problems,  we  advocate  the  use  of  a more  general  spatial  approxi- 
mation method,  an  operator  compact  implicit  method  suggested  by  B.  Swartz  [24]. 
Essentially,  the  same  basic  ideas  are  involved  and  instead  of  setting  up 
spatial  approximations  for  individual  derivative  terms  one  now  poses  the 
difference  approximation  in  terms  of  the  spatial  operator . This  spatial  approxima- 
tion has  been  previously  derived  by  [17];  however \ the  basic  derivation  and 
implementation  there  proceeds  along  lines  different  from  those  taken  here. 

Another  serious  concern  that  one  has  relates  to  the  stability  character- 
istics of  the  overall  method.  If  the  spatial  operator  is  associated  with 
implicit  temporal  schemes,  as  it  might  be  expected,  a variety  of  unconditionally 
stable  schemes  result  for  the  linear  model.  However,  the  cell  Reynolds  number 
stability  characteristics  are  now  somewhat  more  difficult  to  elucidate 
Although  our  analysis  in  section  V is  incomplete,  all  experiments  to  date 
indicate  that  for  the  operator  compact  implicit  (OCT)  approximation,  there  is  a 
wider  range  of  admissible  cell  Reynolds  number  than  for  the  usual  compact  implicit 
methods . 
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In  our  numerical  studies  of  nonlinear  models  we  have  chosen  to  use  two 
different  approaches.  As  a benchmark,  we  have  taken  the  basic  Crank-Nicolson 
routine  solved  bv  simple  successive  approximations.  Our  second  approach  adapts 
a Lees  type  method  [ 12]  which  does  not  require  temporal  iterations  for  a non- 
linear problem.  This  latter  simple  scheme  has  proven  to  be  very  effective  in 
numerical  experiments. 

What  emerges  from  our  investigation  is  that  a promising  class  of  methods 
can  be  developed  around  the  operator  compact  implicit  method.  On  the  basis  of 
our  experiments  an  OCT-Lees  tvpe  scheme  appears  to  be  very  efficient  and 
reliable.  In  the  future  we  hope  to  resolve  questions  concerning  the  treatment 
of  mixed  spatial  derivative  terms  and  to  more  fully  resolve  the  limitations 
associated  with  cell  Reynolds  number  effects. 


I I . BASIC  DIFFERENCE  EQUATIONS 

The  classical  finite  difference  approach  for  solving  two-point  boundary 
value  problems  of  the  form 


(2.1) 


L(u)  = a(x)  u + b(x)  u = f,  xe[0,l] 

XX  X 


with  u(0),  u(l)  given  is  to  separately  substitute  standard  approximations  for  the 
first  and  second  derivatives  in  (2.1)  and  then  solve  the  resulting  system  of 
equations.  Accordingly,  the  centered  second  order  approximation  for  these  terms 
is 


(2.2) 


(2.3) 


5 U.  u, , , - u,  , 
x 1 _ 1+1  J-l  / s , „,,.2\ 

~2h 2h - (ux)j  + °(h  } 


62U. 
x J 

2h 


k-t+1  2U.  + U,  ^ „ 

-±i j 1 - (u^  + 0(h2) 


where  x.  = jh,  j=0, !,•••,. I and  a,  u(x^)  and  h = 1/J  is  the  mesh  size. 


The  resulting  system  of  equations  that  is  derived  upon  substitution  of 
(2.2),  (2.3)  into  (2.1)  is  tridiagonal,  and  hence  easily  solved.  For  the  case 
of  Dirichlet  data,  there  is  no  need  to  create  fictitious  points  (i.e.  to 
extrapolate  information)  in  order  to  implement  the  scheme.  However,  if  higher 
order  accuracy  is  desired,  the  classical  approach  is  to  enlarge  the  basic  mesh 
star,  i.e.  use  more  points  in  the  discretization.  Again,  for  the  centered  type 
of  approximation  fourth  order  accuracy  is  achieved  by  the  following 

(2.4) 


(2.5) 


U 


[• 

[- 


i<2 


6 U.  U. 
-21.1  „ -1-1- 
2h 


8Vi  + 8U1+1  ~ Uj+2 

12h 


(ux)^  + 0(h  ) 


1_  ,2 

12  6x 


2 

6 U. 
x j 


- U1-2  + 16U1_1  - 30U1  + 16Uj+1  - UH2 


12h 


= (uxx)j  + 0(h  ) 
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Bv  substituting  (2.4),  (2.5)  into  (2.1)  a pentadiagonal  system  of  linear 
equations  is  obtained,  and  it  is  necessary  to  use  fictitious  points  near  both 
boundaries . 


A different  fourth  order  approximation  can  be  obtained  by  following  a 
suggestion  of  Kreiss  [14].  The  resulting  representation  is  of  an  implicit 
nature  in  that  there  are  relationships  among  the  function  and  its  derivative  at 
each  of  three  adjacent  mesh  points.  Because  the  method  achieves  the  highest 
order  accuracy  possible  on  the  smallest  star  it  has  been  called  the  compact 
implicit  method.  For  the  derivatives  considered  above,  following  our  notation, 
one  obtains 


(2.6a) 


I + 


T 


_X 

2h 


= (ux>J 


0(h4) 


(2.6b) 


2h 


IT  = 
1 


U.^  - U.  , 

3+1  j~l 
2h 


I + 


x 

6 


(u  ) . + 0(h  ) 
x J 


(ux),+l  + 4(ux)i  + K^-l  4 

X ^ X_J_i  + 0(h4) 


and 


(2.7a) 


2-i 


1 + 


12 


•1  6 2 

~~  U.  = (u  ) . + 0(h4) 
h2  3 xx  j 


(2.7b) 


TV 


Vi  - T 


. + u. , r 6 

1  ^ = I + yif-  (u  K + 0 

2 12  xx  j 


(h4) 


(u  ) + 10(u  ) . + (u  ) . , , 

xx  J+l xx  j xx  j-1  + Q (h4) 


Equations  (2.6)  and  (2.7)  are  derivable  by  either  a Taylor  series  analysis, 
Hermite  polynomial  interpolation  or  by  thinking  of  (2.4)  and  (2.5)  as  Neumann 
series  representations  (up  to  fourth  order)  of  (2.6)  and  (2.7),  respectively. 

As  a reference  for  these  formulas  in  the  case  of  an  uneven  grid,  see  [ 1 ]. 


By  substituting  (2.6)  and  (2.7)  into  (2.1)  it  becomes  apparent  that  in 
general  it  is  not  possible  to  directly  obtain  a tractable  system  of  equations 
terms  of  Uj  alone.  Indeed,  to  solve  the  resulting  system  one  can  define  new 


variables 

diagonal 


F 'v  (u  ) and  S.  ^ 
j xy  j 

system  of  equations 


(u  ) , and  develop  the  following 
xx  j 

approximating  (2.1): 


3x3  block  tri- 


in 
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u4..  - U.  . F + 4F.  + F. 

-iti ill  _ _J±L = n 

U;  2h  6 


(2.8) 


U...  - 2U . + lT  . S,  + 10S.  + S,  , 

(b>  -i+J- i -ti  - -J±l -i ±i  , „ 

h 


12 


(c)  b . F . + a.S.  = f. 

5 J .1  J J 


where  b^  = b(x.)  and  a.  = a(x.)  and  the  above  equations  hold  for  j=l , 2 • • • , J-l . 
Alternatively,  omitting  F.  and  using  only  l'.,  F^  a 2x2  block  tridiagonal  system 
results  from  using  (2.8) (a)  with 


(2.9) 


U . , , - 2U . + U.  , 
J-n  - J Irl 


a[ViF 

1!'vi  1+1  “i 


F , 

j 


+ 


12 


quations  (2.8)  and  (2.9)  require  more  work  to  solve  them  than  the  second 
order  method,  but  generally  the  higher  order  accuracy  permits  one  to  solve  with 
considerably  fewer  points  to  achieve  a comparable  accuracy.  Moreover,  for 
Dirichlet  data,  no  fictitious  points  are  needed.  Boundary  values  (j=0,J)  are 
required  for  Fj  in  (2.9)  and  for  Fj  and  Sj  in  (2.8).  How  these  values  are  obtained 
in  the  time  dependent  case  is  discussed  in  section  ITT. 


At  this  point  we  preview  some  of  the  results  that  will  be  presented  in 
sections  ITT  and  TV  where  a parabolic  problem  with  a spatial  operator  given  by 
(2.1)  is  considered.  There  it  will  be  seen  that  all  the  usual  ways  of  solving 
implicit  systems  incorporating  compact  implicit  schemes  do  not  provide  a 
generally  successful  method  in  the  following  sense.  There  does  not  appear  to  be 
a way  of  achieving  a scalar  tridiagonal  factorization  for  an  unconditionally 

stable  scheme  when  the  compact  implicit  schemes  are  used  for  (2.1).  However,  by 
using  a different  approach  for  the  spatial  operator  these  goals  are  still 
attainable.  Namely,  we  abandon  our  attempts  to  represent  the  separate  derivative 
terms  in  the  spatial  operator  and  adopt  an  approach  which  looks  for  a relation- 
ship on  three  adjacent  points  between  L(u)  and  the  function  u.  The  resulting 
fourth  order  accurate  relationship  may  be  derived  by  a Taylor  series  development 
and  can  he  represented  in  the  following  equations 


(2.10a)  q^(L(U))j+1  + qJ(L(U))j  + q~(L(U))j_1 


hVi  ^ HU1 4 hVi 
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or 


(2.10b)  *-y-  u (x . ) = L(u).  + 0(h4) 

h 1 1 

where  the  operators  Q and  R are  each  tridiagonal  displacement  operators,  namely 


(2.11a) 
(2.11b) 
and  where 


Q l!  = q.  U...  + qT  U.  + qT  U.  i 

i J J+l  MJ  j 1 J-l 


R U.  = rt  U + r°  U.  + r.  U , 
.1  j 1+1  J J J J-l 


+ 2 
q.  = 6a.  a.^  + h(5aj_1  b.  - 2a.  b.^)  - h b.  b.^ 

q j = 4[15aj+i  aj-i  - 4h(ai+i  bi-i " Vi  aj-i}  • h2  Vi  Vi1 


(2.12) 


q.  = 6a.  a..,  - h(5a.,.  b.  - 2a.  b.,.)  - h*-  b.  b 
J J J+l  J+l  J J J+l  J J+l 

rj  = Itqj(2aj+1  + 3h  b)+l3  + qf)(2aj  + hb^  + qj(2aj-l  " hbj-i^ 

rj  = l[ql(2aj+l  + hb j+l 5 + qj(2aj  " bhj)  + qj(2aj-l  ' 3hbi-l)] 

r®  = - (r+  + rT) 

J J J 


These  relationships  were  first  presented  by  Swartz  [24  ] - Equation  (2.10a) 
retains  the  scalar  tridiagonal  feature  of  a second  order  method  while  not 
requiring  additional  fictitious  points  at  the  boundary.  Note,  in  the  case 
where  either  a(x)  or  b(x)  is  identically  zero,  with  the  other  coefficient 
identically  a constant,  the  usual  compact  implicit  schemes  (either  (2.6)  or  (2.7)) 
will  result.  Because  of  these  characteristics  we  have  adopted  the  terminology 
of  referring  to  (2.10)  as  the  operator  compact  implicit  (OCI)  method.  Note  a 
formula  of  structure  similar  to  (2.10)  - (2.12)  is  presented  in  Appendix  A 
for  the  case  of  an  uneven  grid.  In  that  case  the  method  is  third  order  accurate. 

At  least  symbolically,  we  refer  to  the  inverse  of  Q.  The  determination  of 
when  0 can  be  inverted  is  in  general  a difficult  problem.  In  the  case  of  constant 
coefficients  (a(x)  ? a = const,  b(x)  = b = const)  the  invertibility  of  Q on  l. 

hb  ^ 

can  be  fully  analyzed  by  Fourier  analysis  [23].  Defining  R = — as  the  cell 

-1  c a 

Reynolds  number  then  Q exists  for 

(2.13)  R 1 /12  = 3.464 

c 

(The  invertibility  of  Q on  a finite  dimensional  space  is  harder  to  specify. 

For  the  above  case,  a simple  sufficient  condition  guaranteeing  diagonal 
dominance  leads  to  Rc  < (-3  + /24Q) /4  = 3.195.) 


d. 
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The  above  spatial  approximation  can  be  extracted  from  those  schemes  developed 
by  [17].  ilO!  In  the  context  of  approximating  a time  dependent  parabolic  operator. 
Their  resultant  system  of  equations  is  identical  to  an  OCI  approximation  applied 
to  a parabolic  operator  for  some  particular  choice  of  a temporal  scheme.  However, 
the  present  approach  allows  one  to  develop  a variety  of  combinations.  Moreover, 
it  should  be  observed  that  similar  spatial  approximations  have  been  developed 
under  various  names,  in  particular  Collatz  had  sometime  ago  advocated  such  approaches 
which  he  refers  to  as  "mehrstellen"  methods  f 5]. 

In  order  to  apply  any  of  the  methods  presented  here  to  parabolic  equations, 
it  Is.  necessary  to  associate  a time  integration  method.  These  considerations  are 
taken  up  in  the  following  sections. 


Iil.  ALTERNATIVE  TIME  INTEGRATIONS  WITH  COMPACT  IMPLICIT  SPATIAL  DIFFERENCES 

In  this  section  we  consider  time  integration  methods  to  be  used  in  conjunction 
with  compact  implicit  spatial  differencing  [ (2 . 6)- (2 . 9) ] in  the  model  parabolic 
problem 

(3.1a)  u.  = a(x,i)  u + b(x,t)  u 

t xx  x 

(3.1b)  u(0,t)  = c(t);  u(l,t)  = d(t);  u(x,0)  = f(x) 

For  a discussion  of  compact  implicit  methods  applied  to  the  comparable 
second  order  hyperbolic  problem  see  [ 3],  [ 4 ]. 


In  all  that  follows  the  notation  introduced  in  section  II  is  used.  In 
particular , 


(3.2) 


(3.3) 


2-i  ”1 


I + 


I + 


12 


4-  un 

2 i 


2n 


— un 
2h  Y 


where  u”  ^ u(jh,  nAt)  and  At  is  the  time  step. 

Where  the  methods  presented  here  have  appeared  in  a similar  form  in  other  works 
details  are  omitted  and  appropriate  references  are  given. 


Explicit  Methods 

The  first  class  of  time  discretization  methods  to  be  considered  are  explicit 
methods.  Discretizing  (3.1)  in  the  usual  fashion  yields 


(3.4) 


Ilf 1 - U" 

i j , / -n  pn 

At  3j  Sj  + bj  Y 
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!he  right  hand  side  is  evaluated  by  using  the  value  of  l'?  and  (3.2),  (3.3). 
Boundary  conditions  for  s'?  and  F?  must  he  specified.  This  technique  was 
investigated  by  Hirsh  [ 6l,  Rubin  [221,  [23],  and  Adam  [11. 


This  method  is  first  order  accurate  in  time  and  since  it  is  explicit  in  time 
a restrictive  stability  condition  must  be  imposed  (see  e.g.  Hirsh  [23]).  However, 
equations  formed  bv  using  v.3.2)  and  (3.3)  already  are  of  an  implicit  nature. 
Thus,  essentially  no  extra  computational  work  results  if  a second  order  implicit 
temporal  method  is  used  in  order  to  insure  unconditional  stability. 


Implicit  Methods 


Two  second  order  unconditionally  stable  methods  are  considered.  The  first 
ol  these  is  the  usual  Crank-Nicolson  method 


(3.5) 


un+1  - un 
_j 1 

At 


n+1  n+1  , n „n 

a . S . + a . S . 

-)  J J 1 


, n+1  _n+l  , n n 
b . F . + b . F . 

J J J I 


The  second  scheme  to  be  considered  is  adap'ted  from  a : ncond  order  method  presented 
by  T.ees  [12]  which  used  standard  approximations  for  spatial  derivatives.  The 
Lees  approach,  when  implemented  with  compact  implicit  spatial  difference  approxima- 
tions, results  in 


Un+1  - Un_1  an(Sn+1  + Sn  + Sn_1)  bn(Fn+1  + Fn  + F1?'1) 

J J _ 3 1 J J + J J J 2 

2A  t 3 3 


The  Crank-Nicolson  method  would  seem  to  be  more  advantageous  since  it  is 
only  a two-level  scheme.  However,  observe  that  the  coefficients  in  (3.6)  are 
evaluated  at  the  nt"  time  level,  thus  iteration  would  be  unnecessary  even  if  a or 
b were  nonlinear.  Nonlinear  equations  are  discussed  in  greater  detail  in  section  VI. 

Due  to  the  implicit  nature  of  (3.2)  and  (3.3)  one  must  consider  various 
techniques  for  solving  (3.5)  and  (3.6).  Here  we  limit  our  discussion  to  three 
basic  approaches;  predictor-corrector,  block  inversion  and  direct  factorization. 

The  Implementation  of  the  three  methods  is  similar  for  each  of  (3.5)  and  (3.6), 
therefore  details  are  presented  only  for  (3.5). 


1)  Predictor-Corrector  Method.  The  following  predictor-corrector  approximation 
can  be  used  to  solve  (3.5). 


(3.7a) 


U"+1  - un 


n+1  -j-n+1  n n 
. a . S . + a,  S . 

J - 1 j J J 


At 


+ bn  Fn 

j 1 


(3.7b) 


Un+1  - u"  af 1 Sn+1  + an  Sn  bn+1  F^1  + bn  F° 
j j _ j j .1  j J J i j .1 
St 2 + 5 
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Substituting  (3.2)  into  (3.7)  results  in  the  equations 


(3.8a) 


(3.8b) 


+ At  bn  F? 

3 .1 


_ 

2 -1 

-j 

r 

2-1  n 

T 1 /O  n+l 

I - A/2  a , 

j 

2 

6 

X 

"T1- 

T 4-  X 11 

[_T  + I ai 

[I+^]  ^ 

n+1  -=n+l 

j FJ 


where  > = At/h"  and  f']"1’^  in  (3.8b)  is  formed  by  substituting  U?  ^ into  (3.3). 

The  method  may  be  shown  to  be  unconditionally  stable  and  only  requires  the 
solution  of  tridiagonal  matrices,  however,  there  is  one  serious  drawback.  In 
order  to  obtain  a second  order  in  time  accurate  method  it  is  necessary  to  iterate 
(3.8b)  several  times.  Due  to  this  limitation  the  method  is  not  competitive  in 
terms  of  computing  time.  Predictor-corrector  methods  of  this  type  are  examined 
by  Rubin  [23]. 


2)  Block  Methods.  The  block  tridiagonal  methods  fall  into  two  categories; 

3x3  block  and  2x2  block  inversion.  Equations  (3.5)  combined  with  (3.2)  and  (3.3) 
form  the  system.  By  grouping  the  variables  in  vector  format  where  the  unknown 
vector  is 


(3.9) 


one  obtains  the  3x3  block  tridiagonal  equation 


(3.10a) 

(3.10b) 


I + 


T + 


2n 


12 


-1 


x j ,n+l 


2h  j 
-1  , 2 


Fn+1  = 0 


-f-  un+1  - sf 1 = 0 

2 1 j 


(3.10c) 


nn+l  At  n+1  ^.n+1 

“j  -ilj  ?s 


At  n+1  n+1  Itn  At  , n n At  n n 

~2  aj  Sj  ' Uj  + ~ b b + ~ "j  sj 


Methods  of  this  type  were  investigated  by  Hirsh  [ 6 ] . This  method  is  also 
equivalent  to  one  of  the  variants  of  the  Spline  4 methods  of  Rubin  [22]. 
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Alternatively,  substituting  (3.2)  into  (3.5), 
system  with  (3.3),  and  grouping  the  variables  into 


and  completing  the  resulting 
the  unknown  vector 


(3.11) 


results  in  the  2x2  block  tridiagonal  system  investigated  by  Adam  [ 1] 


(3.12a) 

(3.12b) 


In  using  either  of  the  block  methods  it  is  necessary  to  satisfy  extra 

boundary  conditions,  that  is,  a condition  for  F in  the  2x2  block  system  and 

conditions  for  F and  S in  the  3x3  block  system.  Let  us  illustrate  how  this  is 

accomplished  at  the  end  point  X=0.  For  the  3x3  block  method  three  equations 

must  be  obtained  in  order  to  eliminate  F = u I _ and  = u | _.  Combining 

0 x x=0  0 xx  x=0 

equation  (3.10a)  and  (3.10b)  for  j=2  with  (3.10c)  for  j=l  and  the  independent 
Pade  formulas  (see  Hirsh  [ 6 ]) 


(3.13a)  U"+1  - U"+1  + y (f"+1  + Fn+1)  + yy  (s"+1  - Sn+1)  + 0(h5)  = 0 

0 X *—  0 X ^ 0 X 

(3.13b)  U"+1  - u!1+1  + y (F*+1  + 2f!*+1)  - y-  s"+1  + 0(h4)  = 0, 

u l 3 U 1 o 1 

and  are  eliminated. 

Similarly,  for  the  2x2  block  method  two  equations  must  be  obtained  in  order 
to  eliminate  Fq.  Combining  (3.12a)  and  (3.12b)  for  j=l  and  j=2  with  the  Hamming 
formula  (see  e.g.  Ralston  [18]) 


(3.14) 


8U 


11+1  _ 9Un+l  + pn+1 
1 3 


1 fFn+l  + 2Fn+1  - Fn+1) 
h (F0  1 2 


F^  may  be  eliminated. 

The  boundary  condition  for  the  (3x3)  block  method  was  presented  by  Hirsh  [6  ]. 
The  above  boundary  condition  can  be  used  to  retain  fourth  order  accuracy  in 
contrast  to  Adam's  [1]  second  order  boundary  scheme. 


In  order  to  solve  a block  tridiagonal  matrix  the  Thomas  algorithm  [ 8 ] is 
used.  This  algorithm  is  analogous  to  that  of  a tridiagonal  matrix,  with  multipli- 
cation replaced  by  matrix  multiplication  and  division  replaced  by  multiplication 
of  an  inverse  matrix. 
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3)  Factorization  Method.  Bv  directly  substituting  (3.2)  and  (3.3)  into  (3.5) 
one  obtains 


V 


n+l 


tr. 


(3.15) 


2b 


+ _ 
Ah 


>[ri 

,n+1  r r + \ 15  21  5 2 t?n+1  + hn  fT  + 7 * 21  1 « Un 

1 L 6xJ  x q J L 6xJ  XJ 


5 2 U" 

x j 


Upon  examining  (3.15)  it  is  clear  that  there  is  no  way  to  "unravel"  the  implicit 
operators.  This  fact  has  been  observed  by  several  authors  (Ciment-Leventhal  [ 3 ]) 
and  indeed  has  caused  some  to  abandon  entirely  the  compact  implicit  methods  [ 2 ]. 

Following  the  idea  in  f 3 ],  of  completing  the  product,  by  adding  the  second 
order  perturbation  term,  where  is  the  forward  difference  operator. 


(3.16) 


Atlit 

A At 


T 


L-  4 21  1 JL-  |"i  + I 5 2]  1 -2 
12  x J ^2  6 x J 2b 


-*  un 
2h  Uj 


to  (3.15)  the  following  factored  equation  results 

x 

„-i-l 


(3.17) 


_ 1 n+l 

1 - 2 3j 


[l  + i?  v] 

‘"[t  + TZ  v] 


I - f bV 


where 


- r + ~2  a 


> = At/h“  and  p = At/h. 


r[-Kf  ^ 

;M«.r 


T + 7 b 
4 


■j- 


Denote  the  right  hand  side  of  (3.17)  by  G^.then  the  solution  of  (3.17)  can  be 
obtained  by  introducing  an  intermediate  variable  Zj+^  and  splitting  (3.17)  into 


(3.18a) 


(3.18b) 


_ A n+l 
1 ' 


T 8 un+1 

1 ■ A b) 


I + 


j- , r 

12  x J 

6 2]  5 lU 

6 x J x 


7n+l  = rn 

j 'j 


n+l 

j 


7n+l 
'j  ‘ 


This  technique  is  analogous  to  D'Yakonov's  method  [13]  for  two  dimensional 
problems . 

The  algorithm  (3.18)  still  requires  the  solution  of  more  than  one  tridiagonal 
system,  but  no  extra  iterations  are  necessary  as  in  the  predictor-corrector  method. 
However,  a more  fundamental  difficulty  persists  in  this  formulation.  An  inter- 
mediate boundary  condition  for  is  needed.  This  intermediate  condition  implies 

that  either  u or  on  the  boundary  must  he  given.  However  specifying  these  in 
general  (for  example  by  extrapolation)  could  create  an  ill-posed  problem  and 
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generate  instabilities  (numerical  experiments  conducted  by  the  authors  have 
revealed  such  instabilities).  In  certain  problems  however,  e.g.  boundary  layer 
equations,  from  physical  considerations,  extra  conditions  may  apply, then  (3.18) 
may  actually  be  a reliable  and  efficient  method. 


TV.  THE  OPERATOR  COMPACT  TMPETCTT  METHOD 


In  this  section  we  employ  the  OCI  approximation  to  the  spatial  operator 
given  in  (2.1)  with  the  two  time  discretization  methods  considered  in  section  III 
in  order  to  provide  a method  for  solving  the  time  dependent  parabolic  problem  (3.1). 
The  method  is  then  extended  to  two  dimensional  problems. 


The  methods  presented  here  are  unconditionally  stable.  However,  as  with 
most  other  methods  for  this  problem  there  is  a cell  Reynolds  number  condition 
(2.13).  The  discussion  of  stability  will  be  reserved  for  section  V. 

IV. 1 One-Dimensional  Problems 


Consider  the  equation 

(4.1)  u = a(x,t)u  + b(x,t)u  = L(u) 

t xx  x 


where  the  coefficients  in  (4.1)  depend  on  time.  Let  n indicate  the  time  dependence 
in  the  difference  approximation  to  L(u)  at  the  n^  time  level,  i.e.. 


(4.2) 


(LU)n  = 


(0n)-1  Rn 


un 

h2 


The  first  time  discretization  method  considered  here  is  Crank-Nicolson. 


(4.3) 


Un+1  - un 

J .1 

At 


(O^1)-1  Rn+1  un+1 
J 


(0n)  1 Rn  U" 


2h 


which  requires  that  one  solve 


(4.4)  |l  - >(Qn+1r1  Rn+1Ju"+1  = [T  + x(°n)  lRn]Uj  * 

2 

where  X = At/2h  . (Note  well,  for  simplicity  in  the  presentation  of  the  equations 
we  will  be  redefining  X from  time  to  time.) 


Denote 

(4.5) 


the  righthand  side  of  (4.4)  by  0.,  then  (4.4) 
n+ll.,n+l  _ nn+l  „n 


Qn+1  _ XRn+l 
_ — 


j 


j 


can  be  expressed  as 
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Note  the  following  facts  about  (4.5). 

1)  The  matrix  represented  by  Qn+*  - >Rn+^  is  tridiagonal,  thus  very  easily 
solved . 

2)  No  fictitious  points,  or  extra  boundary  conditions  are  needed. 

3)  The  righthand  side  G^  may  be  computed  by  the  simple  recurrence  relation 

(4.6)  Gj  = 2Un  - g"_1 

4)  The  method  is  second  order  accurate  in  time,  fourth  order  accurate  in 
space,  and  unconditionally  stable  (see  section  V). 


The  second  method  to  be  considered  is  derived  from  a Lees  type  scheme  (see 
section  III).  The  Lees  method  combined  with  an  operator  compact  implicit  spatial 
differencing  suggests  the  following  method. 


(4.7) 


uf1  - un_1  (qn)_1  Rn(Un+1  + U?  + u?_1) 
i J J j j 


2At 


3h 


T 


which  requires  the  solution  of 

n+1 


(4 


.8)  I - A(qn)-1  Rn  t/I?+1  = A (Qn)_1  Rn  U"  + 7 + A (Qn)  1 RnJu"  \ 


where  now  * = Multiply  (4.8)  by  Qn  to  obtain 

31i 


(4.9) 


0n  - ARn  I'1}41  = >Rn  u"  + u?-1  + 0°  U?'1 

L J j L .1  1 J i 


Note  the  following  facts  about  (4.9). 


1)  The  matrix  to  be  solved  is  tridiagonal. 

2)  No  fictitious  points  or  extra  boundary  conditions  are  needed. 

3)  The  righthand  side  is  easily  computed. 

4)  The  method  is  second  order  accurate  in  time,  fourth  order  accurate  in 
space,  and  unconditionally  stable  (see  section  V). 

5)  It  is  necessary  to  generate  uj  by  some  other  method  to  begin  the 
computation. 

6)  No  iteration  is  necessary  for  a nonlinear  problem. 

This  last  point  becomes  very  important  in  many  applications,  such  as  the 
boundary-layer  equations. 
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TV.  2 Two-Dimensional  Problems 


R'e  now  turn 

to  the 

consideration  of  the  two 

d imensional 

parabolic  problem 

(4. 

.10) 

u = L 

(u)  + 

L (u)  “ L(u) 

t 

X 

y 

where 

(4. 

, 11a) 

L (u) 

= au 

+ bu 

X 

XX 

X 

(4. 

,11b) 

L (u) 

= cu 

+ du 

V 

yy 

y 

As 

pointed 

out  in 

the  discussion  of  factorization 

methods  in 

section  III 

our  factorization  technique  can  not  be  properlv  adapted  with  the  usual  compact 
implicit  method  to  spatial  operators  with  different  order  terms.  Thus,  the 
discussion  here  is  restricted  to  the  implementation  of  the  OCT  method. 


For  simplicity  (A. 10)  is  solved  on  a rectangular  region  given  by 

{ (x. ,y  ) ; x.  = jh  ; i=0,l,*-*J,  y = kh  ; k=0,l,*«*K} 

■)  k j x k y 

where  boundary  data  is  prescribed  for  all  t for  j=0,J  and  for  k=0,K,  and  initial 
data  is  prescribed  for  t=0.  As  in  [ 3 ] it  is  possible  to  directly  extend  the 
method  developed  here  to  rectangular-like  L-shaped  domains. 


Denote  the  OCT  approximations  to  the  operators  in  (4.11)  by 
1 n U? 

U = (qV1  e" 

li.k  X *h.2 

and 

(4.12b)  |L  111 

L y J i , k 

The  methods  to  be  presented  are  of  the  ADI  (Alternating  Direction  Implicit) 
variety  and  their  derivations  are  similar  to  those  developed  in  [ 3 ] for  the 
treatment  of  the  wave  equation. 

1 rank-Nicolson  Time  Discretization 


T4 . 12a) 


T 

Jx 


As  before  the  first  method  to  be  examined  uses  a Crank-Nicolson  time 
discretization 


(4.13) 


lif1  ~ U"  , (Qn+1)_1  Rn+1  u"+1  + (0n)_1  Rn  U»  . 

.1 , k j , k wx  x j,k  x x 1,k 


TT 


2h 


(Of1)'1  Rf 1 Uf 1 + (O^)-1  Ry  U»iR 
+ 2h  2 

y 
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As  in  the  one  dimensional  case  where  each  of  the  derivatives  was  represented 
separately,  there  is  no  way  to  "unravel"  the  different  inverse  operators  in 
(4.11).  However,  by  adding  to  (4.13)  the  by  now  familiar  second  order  perturba- 
tion cross  term 


(4.14) 


2 

At t_ 

4 2 

At 


mV1  r "I 

(Qn)  1 R n 

X X 

y v 

0 

2 

h 

h 

L x 

V 

j,k 


where  ' is  the  forward  difference  operator, 
the  resulting  enuations  are  easilv  seen  to  assume  the  factored  form: 

(4.15) 


I - 


X (Qn+1)_1  R 

X X 


n+1 
k 


n+1l  [l  - X (of1)'1  Rn+1luf 
xJL  y y y J 3. 

I + V (Qn)_1  R"1  fl  + i ((A'1  R"V  , 

L X 'x  xj  l y y yj  3 .k 


where  A = 
x 


At 

2h 


7t  and  A = 
2 V 


At 


2h 


2 • 


Bv  introducing  the  intermediate  variable  Z?+1  (4.15)  splits  into  two  tridiagonal 

J » K 

systems 


(4.16a) 

(4.16b) 

where 

(4.17) 


I - A (0n+1)_1  Rn+1 

X X X 

I - A (O^1)'1  Rn+1 

y y y 


7n+1 


un+;  = z n+;  , 

J,k  .1  >k 


g"  = [t  + a (0n)  1 Rnl[l  + A (0n)  1 R^u1?  . 

j,k  [ -x  x JL  y y yj  J»k 


Formula  (4.16)  is  analogous  to  an  ADT  type  approximation  solved  with  a D'Yakanov 

splitting.  0n  is  easily  computed  using  previous  values  by  the  following 
-J 


relationship 


^,k  - 2(f ,k  - 7",k + vo-1  u;,k> + c?;i  • 

In  order  to  solve  (4.16a)  boundary  conditions  for  Z^  1 on  the  x = const. 

boundaries  are  needed.  Likewise,  in  order  to  solve  (4.16b)  boundary  conditions 

for  Zn+1  on  phe  y = const,  boundaries  are  needed.  These  intermediate  boundary 
1 »k  , 

conditions  ai\e  obtained  in  the  following  manner: 
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1)  Use  one  sided  differences  to  compute  Z1?^  at  the  four  corner  points. 

n+1  J » * 

Here,  the  fact  that  Z is  a fourth  order  approximation  to 

J » k 

n+l  At.  . . .n+l  . , 

u.  , r(cu  + du  ) . . is  used. 

J , k 2 yy  y j , k 


2) 


n+l 


On  the  x = const,  boundaries  (A. 16b)  is  employed  to  solve  for  Z , 

1 > k 


0n+1  zf?  = 

v 1 ,k 


Qn+1  - X R"+1 

y y y 


n+i 

Ui,k 


,n+l 


3)  Now  that  the  x = const,  boundary  data  for  Z.  , have  been  obtained,  one 

J 9 K- 

can  proceed  with  the  x sweeps  of  the  ADI  scheme  using  (4.16a).  Included 
in  these  sweeps  are  the  y = const,  boundaries.  Thus,  the  Z1}  ^ boundary 
values  necessary  for  the  y sweeps  in  (4.16b)  are  now  fully  available. 

Lees  Time  Discretization 

Finally,  a method  which  is  a generalization  of  the  one  dimensional  OCI-Lees 
scheme  is  examined.  Approximate  (4.10)  by 


(4.19) 


un+;  - ci 

1 >k  .1  , k _ 
2At 


/rin.-l  Rn  /nn.-l  Rn 

(QJ  _JL_  + (Qv)  _2L 

A 9 y *■ 


* u°.k + 


Again,  in  order  to  obtain  a factored  tri-diagonal  method  one  adds  the  second 
order  perturbation  term 

_ nr  „ 1 

5t  U\ 

~At  J’k 


to  obtain 


(Qn)_1  Rx" 
x 2 

Jo0)'1  R^' 
y r i 

h 

h 

L x 

y 

(4.20) 


\i  - * ^v1  Rnl[i  - x (qv1  Rn>n+; 

L x V xJL  y y yj  j .k 

= ft  (Qn)_1  Rn  + X (qV1  Rnlu" 

Lx"x  X y y yj  j .k 

+ fl  + X (oV1  Rnlfl  + x (0n)-1 

1_  x v x xJL  y y yJ  J ,k 
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Denote  the  righthand  side  of  (4.20) 

z"  * and  apply  a D'Yakanov  splitting 
J » k 


bv  0?  , , introduce  an 

j 

to  obtain 


intermediate  value 


(4.21a) 

><- 

X RnlZ 
x x j , k 

-<£ 

rn 

(4.21b) 

[Q°  - 

x Rn]UI?+* 

y y j.k 

- q" 

?U+l 

y 

y 

i »k 

Note  the  following: 


1)  There  does  not  appear  to  be  any  simple  algorithm  for  computing  the 

righthand  side.  However,  upon  multiplying  G by  Qn  (as  in  (4.21a))  it 

^ ^ n 

is  clear  that  only  a backsolve  of  the  tridiagonal  matrix  for  different 
righthand  sides  is  required. 

n+i 

2)  The  intermediate  boundary  condition  for  Z,  , is  obtained  in  the  same 

J ’ K n+l 

manner  as  in  the  Crank-Nicolson  case  once  the  Z^  ^ at  the  four  corner 
points  are  computed. 

3)  As  in  the  one-dimensional  problem  an  extra  plane  of  information  must  be 
generated  to  begin  the  computation  and  no  iteration  is  necessary  for 
nonlinear  problems. 


V • STABILITY  CONSIDERATIONS 

Tn  this  section  we  discuss  two  stabilitv  characteristics  which  enter  into 
the  evaluation  of  the  usefulness  of  difference  schemes  for  parabolic  equations. 

At  the  threshold  one  must  consider  the  Lax-Richtmver  stability  of  the  evolutionarv 
operator  ri9].  More  recently,  it  has  come  to  he  appreciated  that  the  stabilitv 
characteristics  associated  with  the  spatial  operator  should  be  examined  [201, 
f 71,  fill.  The  ability  of  a spatial  difference  scheme  to  resolve  the  spatial 
variation  in  a reeion  of  sharp  gradients  (boundary  laver)  often  gives  rise  to  a 
so  called  cell  Revnolds  number  condition.  Here  we  examine  these  stabilitv  questions 
for  the  compact  implicit  schemes  previously  discussed. 

V.  1.  Temporal  Stability  Analysis 

For  the  case  of  constant  coefficients  one  can  analvze  the  L?  stabilitv 
of  the  difference  scheme  of  interest  by  Fourier  analysis  fl9 1 . Here  the  discussion 

is  limited  to  OCT  schemes.  Substituting  U™  = e*^  into  f4.3)  yields 


(5.1a) 


where 
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24 (cos0-l)+iR  (12-R2) sinfl 

(5.1h)  ?(*)  = 3a — 

30-2R2+(6-R  )cosP+i 3P  sin6 
c c c 

The  term  0(B)  is  associated  with  the  Fourier  transform  of  the  spatial 

2 

operator  alone  [24],  For  stability  is  is  renuired  that  lp„„|  < 1.  Imposing 

C.N 

this  condition  directly  on  (5.1)  vie]ds,  ReP(R)  1 0 as  a necessary  and  sufficient 
condition  for  stability.  This  latter  condition  requires  that 

24(cos°-l)  f 30-2R2+(6-R2)cosR  1 + (12-R2)  3R2sin20  £ 0 . 
c c c c 

Collecting  terms  and  factoring  out  a (eosB-l)  term  yields 

(5.2)  (coso-l) [72n-R4R2+3R4+cos0(144-6OP2+3Ri)l  5 0 

c c c c 

2 

Regrouping,  and  noting  from  (2.13)  that  the  region  of  interest  is  1 12,  yields 

(5.3)  (cosO-l) ri2(12-R2)+12(12-R2cos0)+2B8+(144-72R2+3R4) (cosO+Dl  £ 0 . 

c c c c 

2 

To  see  that  this  inequality  is  always  satisfied  for  I 12,  note 

c 

that  the  term  in  the  left  parentheses  is  1 0 and  the  term  in  the  bracket  is  the 
sum  of  four  terms,  the  first  three  of  which  are  clearly  non-negative . The  last 

2 

term  in  the  bracket  takes  on  a negative  minimum  at  R =12  and  even  when  cosB=l 

c 

this  minimum  is  just  the  negative  of  the  third  term.  This  establishes  that 

o 

f 1 for  R^  f 12,  and  the  unconditional  temporal  stability  of  OCI-CN. 

To  see  that  OCT-Lees  is  similarly  stable,  substitute  U?  = pf  e*2 

1 L 

into  (4.7)  to  obtain  a quadratic  for  Pj , 

(5.4)  p2  + |(K+l)pt  + K = 0 

2 

where  K = p__,  as  in  (5.1)  above  (with  X replaced  by  -r-  X).  Since  the  OCI-CN 
N 5 2 

method  is  unconditionally  stable,  clearly  in  the  range  of  R~  t 12,  |k|  1 1. 

The  stability  of  the  OCT-Lees  method  is  now  contained  in  the  statement  of  the 
following  lemma. 

Lemma . For  the  roots  p.  of  (5.4) 

|PLI  < 1 iff  | K | 1 1. 
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Proof.  First  we  prove  the  lemma  for  the  case  of  equality  in  both  Inequalities. 

i*  , — i<i>/2  , — +lip 

Sav  K = e then  solve  for  P^  directly  as  P^  = P e where  p = e , 

cos  + = -2cos'K  Clear Iv  such  exists  and  thus  |p  I = |p|  = 1.  On  the  other 

II  i <h  / ? 

side,  if  'p  =1,  sav  0^  = e ~ , then  solving  for  K yields 


K = 


l+2e 
1+2  e' 


ii/2 


■U/2 


Thus  <K>  = 1.  This  completes  the  proof  that  |p  | = 1 iff  ( K ) = 1.  To  show  that 

'o'  < 1 iff  1 < 1 examine  the  variation  of  the  roots  pT  with  respect  to  the 

1 1 
unit  circle  as  K varies  from  0 to  +°°.  At  K=0,PT  = 0,  - 4 .both  roots  are  inside 

the  unit  circle.  Now  by  a connectivity  argument,  and  the  fact  that  the  p^  roots 

depend  continuously  on  the  coefficient  K [15],  varying  K such  that  ! K ) <1  then  the 

corresponding  roots  . ^ must  remain  strictlv  inside  the  unit  circle.  Indeed,  if 

some  root  "touched”  the  unit  circle,  i.e.  ' P.  ' = 1,  then  by  our  proof  above 

K = 1.  This  argument  demonstrates  that  for  !k|  < 1,  | p_ I < 1.  Conversely  at 

V = +«  both  p roots  are  outside  the  unit  circle  thus,  again  bv  a connectivity 

argument  the  p^  must  remain  outside  the  unit  circle  for  all  K such  that  | K I > 1. 

Finally  the  stability  of  the  two  dimensional  l.ees-OCI  method  is  established 
using  the  above  Lemma.  Substituting  l,T}  , = pn  e*  ^ ^ ^ into  (4.20)  one  obtains 

1 .k 

(5.5)  P2  + 1(1  - aB)p  - o8  - 0 

where  at  = v . 8 = L l and  where  KP)  and  f (4>)  are  defined  by  (5.1b) 

l->x? (n)  l->yf ( •) 

for  x and  v,  respectively.  Noting  that  a,P  are  each  separately  in  the  form  of 

a 0_„  as  found  above,  one  concludes  from  the  above  lemma,  that  in  the  range 
CN 

■ Rx ' < /l2  , Ir^I  < Si ~2  (i.e.  where  the  cell  Reynolds  number  invert  ibility 
1 c ' - c 

condition  is  satisfied  for  each  spatial  operator)  [ ot  [ - 1,  |R|  * 1.  Now 

identifying  K = - ap  in  (5.5)  clearly  our  above  Lemma  implies  |p|  < 1. 

V . 2 Spatial  Stability 

Experience  with  computations  involving  diffusion  convection  equations  has 
long  shown  that  nonphysical  oscillations  will  appear  in  the  computed  solution 
when  the  spatial  mesh  size  is  not  sufficiently  small  [20),  [7),  [11).  Here  we 
use  the  standard  linear  analysis  to  attempt  to  predict  some  of  the  cell  Reynolds 
number  limitations  associated  with  the  methods  discussed  in  this  paper.  Through- 
out this  subsection,  for  discussion  purposes,  we  will  consider  the  following 
model  "boundary  layer"  problem 
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(5.6) 

au  - bu 

XX  X 

= 0 , a ,b 

positive  constants 

u(0)  = n , 

u(l)  = 1 

where 

in  general  b/a 

is  large. 

Note,  the  solution 

of  (5.6) 

(5.7) 

u (x  ^ ) = Cj 

*x 

+ c0  e3  = 

+ V 

C1  + c2  e , x. 

= j Ax . 

Operator  Compact  Implicit 

The  spatial  stability  analysis  for  this  method  is  quite  straightforward 
and  provides  a practical  guide  for  the  range  of  usefulness  of  the  scheme. 

Assuming  0 ^RU  = 0 is  applied  to  (5.6)  then  one  is  to  consider  the  three  point 

homogeneous  difference  equation 

(5.8)  RU  = 0 

Substituting  a solution  of  the  form  = p ■*  into  (5.8)  (using  (2.11b))  leads  to 
the  following  general  difference  solution 

24+R  (12-R2) 

(5.9)  U.  = c + c^p1:  p = §— 

24-R  (12-R  ) 
c c 

Three  cases  are  possible  for  general  R : 

c 

1.  R < *T7,  p > 1.  The  difference  solution  is  monotone  increasing, 

c 

concave  up  and  properly  approximates  the  true  solution 

2.  /12  < R < 4.207607  (R  value  where  numerator  of  p vanishes), 

c c 

0 < p < 1.  The  difference  solution  is  monotone  increasing  but  concave 
down  and  completely  wrong. 

3.  R > 4.207607,  - 1 < p < 0.  The  difference  solution  is  oscillatory. 

c 

In  summary,  the  spatial  modal  analysis,  of  essentially  the  operator  R 
indicates  that  the  cell  Reynolds  number  Rc  should  be  restricted  to  the  exact 

same  condition  used  for  the  invertibility  of  0,  i.e.  Rc  t /l2\  This  represents 

no  additional  limitation  on  how  one  would  prudently  employ  the  OCT  method. 

Compact  Implicit-Block  Methods 

To  check  the  spatial  stability  of  any  of  the  block  tridiagonal  compact 
implicit  methods  it  is  sufficient  to  consider  any  one  of  them  since  each  method 
(either  the  2x2,  or  the  3x3)  has  the  same  set  of  characteristic  roots.  Thus,  the 
fundamental  modes  of  the  system  can  be  obtained  by  taking  a solution  to  (2.8  a,c) 
of  (5.6)  of  the  form 
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A nontrivial  solution  results  if  the  determinantal  equation 

(5.11)  (y-1)  [ (4-R  )y3+(12-llR  )u2- (12+1 1R  )y-(4+R  ))=0 

c c c c 

holds.  A studv  of  (5.11)  will  at  least  provide  an  indication  of  what  types  of 
nonphvsical  results  are  possible.  However,  there  are  four  roots  (and 
corresponding  arbitrary  constants)  to  contend  with  now.  A proper  analysis 
involves  consideration  of  the  particular  schemes  used  to  approximate  the 
required  derivatives  at  the  boundaries.  Here  we  present  a qualitative  analysis 
of  the  possible  numerical  solutions  of  (5.6)  along  with  some  illustrative 
computational  experiments. 


For  our  model  example  (5.6)  one  would  like  to  obtain  a U.  which  is 


monotone,  or  at  least,  does  not  have  large  oscillatory  modes  which  are 

dominant.  Generally,  this  is  accomplished  by  restricting  R so  that  if 

c 

Re  y < 0 then  [ y : 1 1.  However,  a simple  inspection  of  the  bracketed  cubic  in 
(5.11)  at  y = -1,0,1  reveals  that  such  a condition  can  not  be  found,  since  there 
are  alwavs  three  real  roots  of  (5.11),  y+,  v , y^  such  that 


u+  > 1,  y_  < - 1,  - 1 < yQ  < 0. 

Thus  the  block  tridiagonal  schemes  for  (5.6)  do  not  satisfy  what  has  been 
generally  considered  as  a reasonable  stability  requirement.  Yet  the  schemes 
are  useful  in  practice,  see  section  VI.  The  reason  why  the  oscillatory  modes 
do  not  even  appear  in  some  calculations,  let  alone  dominate  them,  is  tied  to  a 
consideration  of  the  way  the  coefficients  are  determined  by  the  boundary 
conditions. 


A series  of  numerical  experiments  was  made  for  (5.6)  and  qualitatively  we 
can  conclude  the  following.  In  the  range  of  R values  (0  < R < 4//L5  = 1.0328) 

C ^ 

where  u+  1 fy  ! no  dominant  oscillations  occur.  While  in  the  range  - R^  1 
2.14383  corresponding  to  the  y+  range  |y  | f y+  1 e c the  negative  oscillations 
tend  to  affect  more  of  the  region.  For  y+  > e c the  oscillations  are  apparent 
in  most  of  the  region.  Typical  results  are  presented  for  R^  = 1.0,  1.5,  2.0, 

2.4  in  Tables  5.1,  5.2. 


The  case  R = 2.4  is  particularly  interesting  because  here  y = - 2 and  the 
c _ 

y_  term  is  the  dominant  term  in  the  solution  in  the  interior  part  of  region  as 
is  apparent  by  observing  that  the  ratio  of  successive  terms  is  - 2. 


Since  the  circumstances  where  these  spatial  oscillations  will  dominate 
(they  are  always  present  in  general)  is  not  easily  anticipated,  one  should  be 
aware  of  this  potential  problem  for  the  block  compact  implicit  methods. 
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Table  5.1 


Compact  Implicit  (2x2)  Block  Tridiagonal  Solution  of  (5.6) 


1.0,  b/a  = 30. 


R = 1.5,  b/a  = 45. 
c 


1651 3E 
61  3 7 8E 
1 8 130F 

5 1 4 1 4 F 
14  1 34F 
3 8 5 3 3 E 
1 0 4 « 6 E 
2848I1F 
7739SF 
21  01  OF 

6 7 0 8 7 F 
1 5495E 
421  04F 

1 1428E 
31 053F 
84  ? 7 91 

2 2 9 0 3 F 
62155F 
16892E 
48839F 
1 2458F 

3 3 8 0 6 E 
91  886E 
2493  IE 
67769E 
1 8 38hF 
49982E 
1 356  0E 
36864E 
1 0000F 


-12 
-12 
-11 
-1  1 
-10 
-10 
-09 
-09 
-09 
-08 
-08 
-07 
-07 
-06 
-06 
-06 
-06 
-05 
-04 
-04 
-03 
-03 
-03 
-02 
-02 
-0  1 
-0  1 
♦ 00 
♦ 00 
♦ 01 


. 16079P-12 
.59786F-1? 
. 1 7 86  OF - 1 1 
.501 66F- 1 1 
. 1 3794F-1  0 
•3765HF-10 
. 1 0253E-09 
.27885^-09 
•75816F-09 
.2061 1F-08 
. 5602  7F-08 
. 152  10F-07 
.41399F-07 
. 1 1 254F-06 

• 3059QF-06 
. 33 1 53F-06 
.22603P-05 
.61 442F-05 
. 16  702F-04 
.454 O0F-O4 
. 12341F-03 

• 33646F-0  3 
.91 1 8 8 F — 0 3 
.2478HF-02 
.67379F-02 
. 18316F-01 
.4978  1 F-0 1 
. 1 3534F  + 00 
. 3678HF  + 00 
.10000F+01 
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Table  5.2 

Compact  Implicit  (2x2)  Block  Tridlagonal  Solution  of  (5.6) 


I 


R - 2.0 
c 

i,  b/a  “ 60. 

R - 2.4, 
c 

b/a  - 72. 

j 

UJ 

u(x^) 

i 

U1 

u (x  1 ) 

i 

0 . 

0. 

l 

0. 

0. 

? 

.40527  F.  -1  1 

.5594bF-?5 

2 

. 1 0396E-09 

.53927F-30 

3 

-.60778E-1  1 

.48933F-24 

3 

-. 1318hF-09 

.64837F-29 

4 

. 16219E-] 0 

.35239F-23 

4 

. 344  76E-09 

. 7201 0C-28 

5 

| 

-.32360E-10 

.26094F-22 

5 

-.60936E-09 

.79432F-27 

ft 

•73391E-10 

. 192H7F-21 

ft 

. 12990E-08 

.H7565F-26 

7 

15680E-09 

. 14252F-20 

7 

-.251 78E-08 

.96525F-25 

ft 

. 34427E-09 

. 10531F-19 

8 

.51 159F-08 

. 1 0640F-23 

9 

-.74643E-09 

.7781 1F-19 

9 

-. 10152E-07 

. 1 1729F-22 

1 0 

. 1 627  7F-08 

.57495^-18 

1 0 

.20383E-07 

. 12929F-21 

1 1 

- • 354  0 1 E -0  8 

.42484F-17 

1 1 

-.40ftftftF-()7 

. 14252F-20 

1? 

.77  089E-08 

.31 39 1 F- 1 ft 

12 

.81 453F-0  7 

. 1571 OF-19 

13 

-. lft777E-07 

.231 95F -15 

1 3 

-.  16283F.-06 

. 1731 7F-18 

14 

• 36522E - 0 7 

. 17139F-14 

14 

. 32573E-06 

.19089F-17 

15 

-. 79496E-07 

. 12664F-13 

15 

-.651 38E-06 

.21 042F-16 

1ft 

. 1 7304E-06 

.9357ftF-13 

1ft 

. 1 3028E-05 

.23195E-15 

17 

-•37ftft7E-06 

.59144^-12 

1 7 

- . 26056E-05 

• 25569F- 1 4 

1ft 

.81 99 1 E - 06 

.51091F-1  1 

1ft 

.521 13E-05 

. 2R  1 H5F-  1 3 

19 

-. 17A47F-05 

.37751F-10 

19 

-. 10423E-04 

.31068F-12 

20 

. 3Aflft 1 E-05 

.27895F-09 

20 

.20845E-04 

• 34  24  7F- 1 1 

21 

- . 8454 1 F- 05 

• 20ft 1 2F-08 

21 

-.41 690E-04 

•37751F-10 

2? 

. 18423E-04 

. 15230F-07 

2? 

.833R1E-04 

.4161 4F-09 

23 

-. 39949E-04 

. 1 l?54F-0ft 

23 

-. 16676E-03 

• 4 5fl  7 2F- 08 

24 

. 88Q83E-04 

1 .83153 F -06 

24 

.33357E-03 

,505bpF-07 

25 

-. 18344E-03 

• ft  1 442F-05 

25 

- . 6665 1 E- 0 3 

•55739F-06 

2ft 

.4ft035E-03 

• 454  0 OF -04 

26 

. 1340  IE-02 

.61442F-05 

27 

-.55?5ftE-03 

• 33546F-03 

27 

-.26014E-02 

.67729F-04 

2ft 

.45I28E-02 

,247«8F-02 

28 

.60827E-02 

. 74659F-03 

29 

. 14551F-01 

. 18316F-01 

29 

-•23290F-02 

.R2297F-02 

. 14782E+00 

. 1 35  34F ♦ 0 0 

30 

. 1 1462E*00 

.9071BF-01 

31  1 

. 1 0000E* 01 

. 1 0000F.01 

3) 

. 1 OOOOE*Ol 

. 1 0000r»01 

! 


24 


7 / - 2 9 


' 


j 

i 

i 

< 

V 

I 


VI.  NUMERICAI,  EXPERIMENTS 
VI. I Introduction 

In  this  section  results  of  numerical  experiments  conducted  with  the 
various  schemes  that  were  discussed  in  sections  IT  - IV  are  presented. 

These  calculations  were  performed  in  order  to  determine  the  viability  of  the 
OCI  method  for  solving  viscous  flow  problems  and  to  understand  its  character- 
istics and  limitations  and  to  compare  its  performance  with  the  classical 
second  order  techniques  nox^  in  general  use  as  well  as  to  other  fourth  order 
approaches . 

One  of  our  major  concerns  is  the  efficiency  of  the  various  schemes,  i.e. 
computation  time  required  to  obtain  a given  accuracy.  Obviously  this  is 
machine  as  well  as  programmer  dependent.  In  order  not  to  bias  any  of  the 
techniques  care  was  taken  to  program  the  algorithms  in  an  efficient  and 
consistent  manner.  The  computing  times  that  are  given  include  time  for: 
matrix  setups,  inversions,  boundary  condition  evaluations  and  (for  nonlinear 
problems)  iteration  procedures.  All  results  were  computed  on  the  NSWC/WOL 
CDC  6500  computer. 

The  operation  count  estimates  (multiplications  and  divisions)  for  the 
block  tridiagonal  inversion  algorithm  is  given  in  [ 8 ] as 

(6.1)  ops  = (3n-2)(n3  + m^) 

where  m is  the  order  of  the  block  and  n Is  the  number  of  equations.  This 
estimate  assumes  full  blocks.  However,  if  the  specific  values  of  the  elements 
of  the  blocks  are  cafcen  into  account,  e.g.  zeroes  and  ones,  the  actual  opera- 
tion count  can  be  greatly  reduced.  A comparison  of  operation  counts  for  the 
various  inversion  procedures  (assuming  full  blocks)  and  the  modified  algorithms 
are  presented  in  Table  6.1.  Also  included  there  are  the  counts  for  matrix 
setup  operations.  Note  that  for  the  block  methods  the  inversion  of  the  matrix 
is  the  dominant  factor  in  the  running  time,  while  for  the  OCI  technique  the 
matrix  setup  accounts  for  most  of  the  time. 

VI. 2 Linear  Parabolic  Equation 

The  first  numerical  experiment  involved  the  solution  of  a one  dimensional 
linear  parabolic  partial  differential  equation  with  variable  coefficients 

(6.2a)  ut  = a(x,t)uxx  + b(x,t)ux  ; t t 0 ; 05x11, 
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where 


a (x , t ) 


1 (x+1) 

2 2 ’ 
(t+2  ) 


b (x , t ) 


I W) 

ft + 2) 


with  the  exact  solution 


(6.2b)  u(x,t)“u  (x,t)»  exp f (x+1 ) (t+2) ] 

Initial  and  boundary  conditions  are  Riven  by 

u(x,0)  = u (x,0) 

(6.2c)  e 

u ( 0 , t ) = u (0,t)  ; u ( 1 , t ) = u ( 1 , t ) 

e e 

Ibis  example  was  constructed  in  order  to  test  the  stability  and  convergence 
properties  of  the  methods  under  consideration  for  a variable  coefficient  problem. 
Results  are  shown  in  Table  6.2  and  Figure  1.  All  the  methods  tested  were 
stable  and  show  the  predicted  convergence  rates.  Crank  Nicolson  temporal 
integration  was  used  for  all  the  schemes. 

Of  basic  interest  is  the  savings  that  can  be  obtained  in  storage  and 
computational  time.  As  noted  in  Table  6.2  and  Figure  1 the  OCI  technique  is 
the  most  efficient  scheme  of  the  methods  tested.  This  is  not  wholly  unexpected 
since  the  block  methods  require  additional  work  to  compute  the  first  and/or 
second  derivatives. 


It  is  also  important  to  note  the  differences  in  the  computed  L9  errors 

of  the  fourth  order  methods.  These  result  from  several  factors  among  which  are 
the  local  truncation  error  and  boundary  conditions.  The  spatial  truncation 
errors,  which  are  dominant  for  the  case  considered,  are  given  below. 

Compact  Implicit  - (Block  Methods) 


First  derivative 


Ef  = 


- h 
180 


0(h6) 


Second  derivative 


Es  = "Z 


h4  vi  n,.6. 
40~  u + °(h  > 


Thus  for  equation  (2.1)  with  a and  b constant  the  local  spatial  truncation  error 
at  point  j would  be 


(6.3) 


. 4 r a v i 
= ~ 6 Uao  u 


240  U1  180 


VT 

UJ 
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OCI 


Specializing  Fq.  (A. 16)  for  constant  coefficients  yields 


(6.4) 


r a vi  , 3b  v\ 

'200  uj  + 200  ’V 


In  achieving  a scalar  tridiagonal  system,  the  OCI  technique  leads  to  an 
unsymmetric  difference  formula  and  thus  has  a larger  local  truncation  error  than 
the  block  methods  that  were  derived  from  symmetric  formulations.  Were  it  not 
for  the  different  boundary  conditions,  Pade  relations  (3.13)  for  the  3x3  blocks 
and  a Hamming  type  formula  for  the  2x2  blocks  (3.14)  both  block  methods 
would  give  identical  errors. 


''T.2.1  General  Boundary  Conditions 

The  OCI  method  can  also  be  applied  to  problems  with  more  general  boundary 
conditions  of  the  form 

(6.5)  A u^  + Bu  = g 

A linear  fourth  order  accurate  expression  is  sought  relating  u^  at  the 
boundary  with  u and  L(u)  at  points  .1*0, 1,2,  i.e. 


16.6)  Fq  = (U„)0  ■ HqUq  + BjOj  + H2Uj  + c0f0  + 

F.mploying  the  differential  equation 


L(u)j  = aS  + bF  = f (j=0, 1,2)  , 


and  the  compact  implicit  formulas 


F0  + 4F1  + F2  ■ h <"2  - U0> 

S0  + 10Sl  + S2  - i|  (U0  - 2Ul  + U2) 

n 


s0  + 4S1  + s2  = h (F2  ~ F0)  ’ 


the  coefficients  in  (6.6)  can  be  evaluated.  These  coefficients,  the  truncation 
error,  and  the  extension  to  time  dependent  problems  are  given  in  Appendix  B. 

As  an  example,  Eq.  (6.2a)  was  solved  with  the  following  boundary  conditions 


0 x = 0 , 

0 x - 1 , 


u + u = u (0)  + u (0)  = (t+3)  exp  [ t+2] 
x e e 

x 

u = u (1)  = exp  [2 (t+2)] 
e 
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Table  6.3  shows  the  I.,,  errors  and  L?  rates  for  different  mesh  widths. 
Comparisons  with  the  results  in  table  6.2  Indicate  that  for  general  boundary 
conditions  the  L,,  error  Is  larger  and  the  computation  time  is  increased. 


VI. 3 Burgers  Equation 

The  results  just  presented  although  extremely  promising  are  obtained  for 
a linear  equation.  In  order  to  investigate  the  various  methods  for  a nonlinear 
problem  that  might  be  indicative  of  viscous  flows  the  one  dimensional  Burgers 
equation  was  chosen.  Consider 

(6.7)  u = - (u-oi)u  + vu 

t x XX 

With  the  exact  steady  state  solution  given  by 


(6.8)  u^(x)  = ajl  - tanh  (^)| 

Near  x = 0,  u(x)  exhibits  large  gradients,  and  as  v -*  0 a steep  shock  wave  forms 
The  ability  to  resolve  this  flow  field  would  demonstrate  the  viability  of  the 
various  methods. 

Solutions  were  obtained  in  the  domain  - 5 f x 1 5 with  a = 1/2  and  for 
various  values  of  v,  and  with  the  exact  values  of  u(x)  specified  at  the 
boundaries.  The  initial  conditions  employed  for  all  cases  are 


u(x,0) 


1 - 5 < x < 0 

.5  x = 0 

0 0 < x < 5 


Results  of  computations  with  the  OCI  (Crank  Nicolson  and  Lees)  methods 
and  the  second  order  Crank  Nicolson  finite  difference  scheme  are  presented  in 
Tables  6.4,  6.5,  6.6  and  6.7  and  Figure  2. 


Since  Eq.  (6.5)  is  nonlinear,  iteration  is  necessary  for  the  Crank-Nicolson 

temporal  discretization.  We  adapt  the  OCI  method  with  successive  approximation 

for  the  nonlinear  term,  tiu  , i.e. 

x 

(6.9)  U"+1(Ux)1+1  ' Uj(Vl+1 

where  Uj  is  the  latest  iterant  value.  This  procedure  converges  linearly. 

The  second  order  finite  difference  scheme  uses  a different  type  of  lineariza 
tlon,  i.e. 
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(Tn+  - n)  {(!n+1  - Un+1  irn  - If" 

(6.10a)  (U  - a)n+  (U  )n+'  = 2- { Jn'2V  J~~  + -•1~1 

1 x j 21  2 Ax  2 Ax 


where  U?+  Is  replaced  by 


(6.10b) 

k 

lr 4 being  the  lastest  iterant, 
convergence  properties  [16]. 


(U*  + u")/2  , 


This  form  of  iteration  has  super-linear 


n+l  * 

Both  methods  assume  an  initial  guess  for  U.  = U.  which  is  used  to 

J j 

solve  the  resultant  tridiagonal  system  of  equations.  Iteration  is  employed 
until  the  difference  between  successive  iterants  is  less  than  some  preset 
tolerance.  The  steady  state  is  assumed  when  differences  in  solution  values 
at  two  time  steps  is  less  than  some  predetermined  value. 

In  contrast  to  the  above  procedure,  the  OCI-Lees  discretization  does  not 
require  iteration  and  generally  approached  the  steady  state  in  about  the 
same  number  of  time  steps  as  the  OCI-CN  method. 

Figure  2 presents  a graph  of  the  computed  error  versus  the  number  of 

intervals,  for  the  fourth  order  and  second  order  schemes.  The  storage 
savings  possible  with  the  OCT  method  are  readily  evident  from  the  figure. 

Tables  6.6  and  6.7  compare  solution  values  obtained  from  the  4th  order  and 
second  order  methods  with  the  exact  value,  for  two  cases,  v = .5  and  v = .031. 

Although  the  cell  Reynolds  number  analysis  for  the  OCI  method 
given  in  section  V was  derived  for  a linear  spatial  operator,  this  theory 
can  be  useful  to  predict  the  behavior  for  nonlinear  time  dependent  problems. 

For  the  Burgers  equation  it  was  found  that  physical  solutions  were  obtained  for 
a steady  state  only  when  R 1 < 2.55,  where 


(u-a) Ax  _ Ax 

v 2v 


However,  a careful  inspection  of  the  numerical  results  indicates  that  for 
fR  [ > 2.55, in  computing  the  transient  solution,  values  are  obtained 

max 

which  yield  cell  Reynolds  numbers  exceeding  /T2 , and  physical  steady  state 
solutions  can  not  be  obtained.  These  results  suggest  that  when  the 
homogeneous  case  maintains,  one  monitor  the  evolution  of  the  local  cell 
Reynolds  number  and  consider  modifying  the  spatial  mesh  when  necessary. 


i 

I 
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In  contrast  to  the  above  behavior,  for  the  boundary  laver  equations 
oscillations  when  they  occurred  were  confined  to  some  local  region,  but 
phvsical  solutions  were  still  obtainable  elsewhere  (see  section  VT.5). 

The  results  of  the  computations  presented  above  indicate  that  the  OCI 
method  can  be  adapted  to  handle  nonlinearities  with  very  little  additional 
effort  and  can  resol've  regions  with  sharp  gradients. 

VI. 4 Two  Dimensional  Problems 

The  OCI  method  was  tested  for  a two  dimensional  parabolic  equation 
(6.11a)  u{  = a(x,y,t)uxx  + b(x,y,t)u^  + c(x,y,t)uv^  + d(x,y,t)u^ 

whose  coefficients  were  constructed  in  order  to  obtain  the  solution 
(6 . lib)  u(x,v,t)  = exp( (x+1) (y+1) (t+l) } . 

Neither  efficiency  studies  nor  comparisons  with  other  methods  were  made. 

The  aim  here  was  mainly  to  check  the  order  of  accuracy  and  the  viability  of 
the  ADI  formulation.  Table  6.9  demonstrates  that  the  splitting  technique 
given  in  section  IV  yields  fourth  order  accuracy.  Ciment  and  Leventhal  [ 3 ] 
have  demonstrated  that  for  hvperbolic  equations  this  type  of  ADI  scheme 
retains  fourth  order  accuracy  on  other  than  rectangular  domains,  e.g. 

L shaped  domains.  Similar  results  are  expected  for  parabolic  equations. 

VI.  5 Boundary  Laver  Equations 

The  main  thrust  of  this  work  is  to  develop  methods  that  could  efficiently 
solve  viscous  flow  problems.  Although  the  next  two  examples  are  rather 
idealized,  they  do  possess  the  intrinsic  features  of  more  complicated  boundary 
layer  flows.  The  two-dimensional  laminar  incompressible  boundary  layer  along 
a flat  plate  with  and  without  pressure  gradient  require  that  one  solve 

(6.12)  uu  + vu  = u u +vu 

x v e e yy 

x 7 7 

u + v = 0 
x y 

with  boundary  conditions 

u(x,0)  = v(x,0)  = 0 ; u(x,°°)  = u , 

e 

and  initial  conditions 

u(0,y)  ■ F(y) 

By  making  the  following  transformation 

■ x , q = y/ /2vx  , v = v/2£/v 
<J>  “ qu  - V 
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the  governing  equations  reduce  to 

(6.13a)  2f,uu  = $u  + u + 2f,u  u 

q no  e ef 

(6.13b)  $ = 2?,u  + u 

n t 


For  the  case  of  zero  pressure  gradient,  ug  = 1,  the  Blasius  flow  is  obtained 

and  the  equations  become  independent  of  r . The  solution  of  the  ordinary 
differential  equation  that  is  recovered  is  compared  with  the  time  asymptotic 
solution  of  equations  (6.13a)  and  (6.13b). 

Due  to  the  nonlinear  term  2tuuf.  and  the  decoupling  of  the  momentum  and 

continuity  equations  iteration  is  required  for  the  Crank  Nicolson  scheme. 

The  term  2,ruu_  is  discretized  as  follows 


un+l  + rn 

(2U'.U.r)n+  = 2^n  ( -1— - 


un+1  - un 


-1  )•(  J — 77 — *■  ) -V  K1  1 


F.mploving  quasi-linearization  (Newton-Raphson  iteration) 

(u^1)2-  2F*  r"+1  - (u*) 2 , 

•k 

where  U , is  the  latest  iterant,  the  following  relationship  is  obtained 


(2nyijr) 


n+* 


& r«.;  vf1  - K)2  - K)21 


i j 


Care  must  be  taken  in  the  choice  of  the  initial  U . Although  setting  U.  = U. 

n+i  2 ' J J 

gives  a second  order  approximation  for  (u  . )“,  it  reduces  to  first  order  for 

• n-t-1  ' 

(2r,uuf)  due  to  the  A6  in  the  denominator  of  the  derivative  approximation 

for  u,. 


Tf,  however,  U is  approximated  by  the  extrapolation  formula 


U*  = 2Uj  - Uj  1 


,n+’ 


second  order  accuracy  is  recovered  for  the  term  (2F,U^U^)  . This  procedure 

was  employed  in  the  Crank  Nicolson  calculations  for  both  second  and  fourth 
order  methods. 
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Since  the  accuracy  of  this  1 inearization  is  of  the  same  order  as  the 
temporal  discretization  error  further  iteration  should  be  unnecessary. 
However,  the  decoupling  of  the  continuity  and  momentum  equations  require 
iteration  to  obtain  the  desired  temporal  accuracy. 

If  4>u  + u is  identified  as  L(u)  and 

n nn 

f(Qn)_1  R^lu1  -v  L0!!1  = A1  + U1 

L J nn  n 


is  either  the  fourth  order  OCI  approximation  or  the  second  order  centered 
difference  approximation  (note  that  in  this  case  Q f T)  and  i and  m denote  appro- 
priate time  levels,  then  the  boundary  layer  equations  may  be  discretized  in  two  ways. 


lie  t hod  I. 
(6.14) 


<vr 


(0n+1  )_1  -^+ 

|un+1  + un' 

.1  .1 

2 

An"  . 

where  g(f,)  is  the  pressure  gradient  term.  Employing  Newton  Raphson  iteration 
and  performing  the  indicated  operations,  the  following  system  of  equations 
is  obtained 


(6.15) 


Qn'H(2?U*)  - \ Rn+  rn+1  = nn+'5  (run)  ^ Rn+  un 

j 2 J 3 |_  J 2 J J 


+ 0n+1<  fu*  + Afg  «>] 


+ AEgCOJ  j=l,-»*,J-l 

where  X = A£/An^  and  £ is  evaluated  at  n+S . The  system  is  tridiagonal  and  can 
be  inverted  given  U at  j=0  and  J. 

Since  <f> j appearing  in  0n+!:  and  Rn  p must  be  evaluated  at  rtb1-  the  continuity 
equation  becomes 

(6.16)  25(1^)"^  + uj**5  - (*n)"+V 


or 


, ^ (of1  - of  (of1  + o")  . 

<‘n>r  ’ 25"  + J'  2 'V°  ' 


which  can  be  integrated  using  Simpsons  rule, 

(6.17) 


<1  ■ + h [Vi  + iV> + vi-i]n+1  • 
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U"'  ■ (Vr)f  = U_  W+V1  kn+1]u"+1  +[(nn)-1Rn]u">  + g (rn+ 

2An  ’ J ' " 


Performing  the  indicated  operations,  the  following  tridiagonal  system  of 
equations  is  obtained. 


(6.19) 

where 


|qn+1(2rn+'  U*)  - jRn+1Ju^ 


n+l  _ ^n+1 


[o"  + cn+:  r*  + atetc"'1*5)] 


c"  - |[(Qn)-1R"]l."  + e"*1  u”  o"  . 


0^+  can  be  evaluated  at  the  new  time  level  by  the  recursive  relations 


(6.20)  GI?+1  = (rti+3/2  + 2fn+1>)  (u"+l)2  - g"  - AFg(Fn+‘  ) - F^d'")2 


The  handling  of  the  ur  term  in  the  continuity  equation  for  the  above 

method  deserves  special  attention.  Since  <J>  and  u are  evaluated  at  the  (n+l) 
time  level,  u,  must  also  be  evaluated  there.  It  was  thus  necessary  to  use 

a one  sided  second  order  accurate  derivative  approximation  for  ur  (uniform 
mesh) 

(6.21)  (ur)i+1  = (3,J1?+1  - 4U"  + U?_1)/2AF. 

J J J J 

Both  methods  I and  IT  worked  successfully.  However,  method  I was 
preferred  computationally. 

Lees  Method 


The  Lees  method,  which  does  not  require  iteration,  was  also  used  to 
obtain  solutions  of  the  boundary  layer  equations.  The  discretized  equation  is 


2 F U 

j 


(6.22) 

which  reduces  to 


„,n+l  „n-l.  , 

n „n  (UJ  - Vi  } = (Qn)~'~  Rn 


2A  F 


3An 


[uj  1 + u”  + Uj+1]  + g(Fn)  , 


(6.23) 


[o"(tn  uj)  - [»”  + U^1] 

+ on[rn  o"  i'"_1  + 6tg«")J 


where  X = AF/An  . 
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In  order  to  update  Q and  R at  the  (n+1)  *■  £ level,  <1  must  first  be 
evaluated  at  that  level.  Hence  it  Is  necessary  to  employ  (6.21)  for  the 
integration  of  the  continuity  equation. 

The  Lees  method  incorporates  three  levels  of  information,  therefore 
a starting  technique  such  as  Crank  Nicolson  is  needed.  Furthermore,  to 
employ  variable  A£  requires  a restarting  procedure  so  that  for  such  cases 
the  Lees  method  may  lose  some  of  its  appeal. 

Table  9 compares  the  computed  L2  error  and  wall  shear,  t , for  the  various 
methods.  Note  that  using  formula  (6.6),  especially  for  the  cases  with 
few  mesh  points  yields  better  predictions  for  the  wall  shear  than  the 
standard  fourth  order  one  sided  difference  formula.  Figure  3 which  is  a 
plot  of  the  L2  error  versus  the  number  of  intervals  illustrates  the 
savings  in  storage  that  can  be  obtained  by  the  fourth  order  methods. 

Effects  of  cell  Reynolds  number  were  noted  for  the  OCI  and  2x2  block 
methods  for  the  case  Aq  = 1.  In  the  discussion  below,  reference  will  be 
made  only  to  OCI  since  the  cell  Reynolds  number  characteristics  of  the 
block  methods  are  not  easily  understood. 

The  local  cell  Reynolds  number  for  the  Blasius  problem  is  defined  as 

R = (JiAn  . 
c 

The  function  $ grows  linearly  for  large  n,  so  that  at  some  point,  as  the 
domain  is  extended,  the  local  value  of  Rc  will  exceed  /T?  and  when  A. 207 
is  exceeded,  as  predicted  by  the  linear  analysis,  oscillatory  behavior  will 
ultimately  result.  This  conclusion  indicates  ironically,  that  in  the 
region  of  small  gradients  (at  the  boundary  layer  edge)  increased  local 
resolution  might  be  required  to  remove  the  oscillatory  behavior.  Keller  [ 9 ] 
has  made  similar  observations  for  second  order  methods. 

In  the  above  calculations  for  Aq  = 1 , with  the  boundary  layer  edge 

at  q = 6 no  oscillations  occurred  (with  the  internal  (R  ) = 3.790  at 

c max 

q = 5.0).  However,  with  the  boundary  layer  edge  extended  to  q = 10.0, 
oscillations  appeared  from  the  outer  boundary  inward  to  q = 8.  Since 
the  velocity  oscillates  (<  1.5%)  about  a value  greater  than  1,  in  the 
region  A 1 q < 10,  the  velocity  profile  exhibits  overshoots. 

Enlarging  the  domain  does  not  cause  the  oscillations  to  invade  the 
region  where  the  cell  Reynolds  number  condition  is  satisfied  as  inspection 
of  the  solution  reveals.  Furthermore,  the  computed  wall  shear  is  affected 
only  slightly  (see  below). 


T ( Eq . 6.6) 

w 

6.0 

. A771 

10.0 

. A869 

20.0 

. A893 
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The  second  boundary  laver  example  considered  was  the  flow  with  an 
adverse  pressure  gradient,  the  Howarth  problem,  where  the  external  velocity 
is  linearly  retarded 

u = 1 - C . 
e 

As  a consequence  of  the  adverse  pressure  gradient,  separation  (vanishing 
of  the  wall  shear)  will  occur  at  some  point  downstream  along  the  flat 
plate.  The  determination  of  the  entire  flow  development,  from  the  initial 
Blasius  profile  to  the  point  of  separation  is  sought.  Since  the  boundary 
laver  equations  break  down  at  separation,  it  is  expected  that  the  numerical 
computations  will  be  sensitive  near  that  point.  Factors  such  as  AC  step 
size  and  iteration  criteria  may  strongly  influence  the  calculations.  A 
complete  analysis  of  the  behavior  of  the  numerical  solution  near  the 
separation  point  is  not  considered  in  this  report.  However,  a set  of 
calculations  for  a fixed  AC  = 10-^  and  iteration  convergence  parameter, 
e = 10" ^ have  been  obtained  and  are  shown  in  Table  10.  The  computed 
separation  point  (the  point  where  the  shear  changed  sign)  along  with  the 
running  time  for  each  calculation  is  given. 

In  general,  one  iteration  was  required  for  convergence  of  the 
Crank  Nicolson  discretization,  except  near  separation.  Several  calculations 
did  not  converge  near  the  separation  point  and  were  thus  terminated  there. 
These  are  indicated  by  asterisks  in  Table  10. 


Hirsh  [ 6 ] employing  the  3x3  block  method  computed  the  separation 

point  at  r = .119818  for  An  = .2.  Since  no  discussion  of  the  behavior  of  the 

velocity  profiles  near  separation  was  given  in  [ 6 ] a detailed  comparison 

can  not  be  made.  However,  our  calculations  show  that  the  second  order 

Crank  Nicolson  and  the  2x2  block  methods  give  £ > .1200  for  all  An, 

sep 


whereas  OCI  (in  particular  the  Lees  discretization)  yields  values  of  ( < .1200. 

sep 


Near  separation  the  cell  Reynolds  number  becomes  very  large  and  exceeds 
the  limits  set  by  the  linear  theory.  However,  the  actual  behavior  of  the 
solution  in  this  region  is  not  predicted  well,  and  simple  monitoring  of 
Rc  is  not  helpful  as  it  was  for  the  Blasius  problem  since  it  is  now  necessary 
to  take  into  account  the  nonhomogeneous  terms  contribution.  Further  investiga- 
tion is  necessary  to  understand  even  the  nonhomogeneous  linear  case. 
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Table  6.1 

Matrix  Setup  and  Inversion  Operations*  Uniform  Mesh 


MATRIX 

INVERSION 

SETUP 

TOTAL 

ESTIMATED** 

ACTUAL 

Scalar 

Tridiagonal 

(OCT-CN) 

5N-4 

5N-4 

22N-22 

27N-26 

2x2  Block 

Tridiagonal 

(C-N) 

36N-24 

27N-60 

8N+16 

35N-44 

3x3  Block 

Tridiagonal 

(C-N) 

108N-72 

49N-62 

4N+24 

53N-38 

Scalar 

Pentadiagonal 

11N-16 

ION*** 

21N-16 

* Here  it  is  assumed  that  multiplications  and  divisions 
are  equivalent.  However  on  certain  machines  this  may  not  be 
true,  e.g.  on  the  CDC  6600  a division  is  comparable  to  six 
multiplications.  The  operation  counts  would  have  to  be  changed 
accordingly  for  the  methods. 

**  Reference  [ 8 ] ■ 

***  Does  not  include  extrapolation  formulas  for  points  adjacent 
to  the  boundaries. 


36 


NSWC/WOL/TR  77-29 


Table  6.2 

Linear  Variable  Coefficient  Parabolic  Equation 


u = a(x,t)u  + b(x,t)u 
t xx  x 

u = exp{ (x+1 ) (t+2) } 


TIME 

COMPUTING  TIME* 

METHOD 

N 

STEPS 

L2 

ERROR 

l9  rate 

(SEC) 

Second 

100 

2000 

.20 

* IQ’04 

35.5 

Order 

Crank 

160 

2000 

.79 

* io-°5 

1.93 

55.4 

Nicolson 

200 

2000 

.51 

* 10-05 

1.96 

68.8 

* 10-05 

1.97 

400 

2000 

.13 

135.5 

3x3 

5 

2000 

.15 

* 10-°4 

12.2 

Block 

-06 

4.06 

17.7 

Crank 

10 

2000 

.90 

* 10 

Nicolson 

20 

2000 

.51 

* 10-07 

ln~08 
* 10 

4.14 

4.67 

30.4 

40 

2000 

.20 

58.5 

-05 

2x2 

Block 

5 

2000 

.83 

* 10 
* 10-06 

3.58 

7.4 

Crank 

10 

2000 

.70 

11.9 

Nicolson 

20 

2000 

.48 

* 10-07 

in-08 
* 10 

3.87 

4.59 

20.8 

40 

2000 

.20 

40.2 



-04 

Operator 

5 

2000 

.24 

*10 

5 . 4 

Compact 

-05 

4.00 

Implicit 

10 

2000 

.15 

* 10 

8.8 

Crank 

Nicolson 

20 

2000 

.94 

* 10-07 
* 10-08 

4.00 

4.52 

15.6 

40 

2000 

.41 

30.0 

* Computation  times  are  for  a CDC  6500 
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Table  6.3 

Linear  Variable  Coefficient  Parabolic  Equation 

u = a(x,t)u  + b(x,t)u 
t xx  x 

u = u = exp(x+l)  (t+2) 
e 


u (0)  + u (0)  = (t+1)  exp  [t+2]  , u(l)  = u (1) 
x e 

OCT  - 2000  TIME  STEPS 


N 

l2  ERROR 

l2  RATE 

COMPUTING  TIME* 
(SEC) 

5 

.222  * 10"02 

6.57 

10 

.122  * 10"03 

4.796 

9.58 

20 

.737  * 10-05 

4.049 

16.71 

40 

.441  * 10-06 

4.063 

30.42 

* Computation  times  for  a CDC  6500 


i 

i 


■J 
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Table  6.4 

Steady  State  Solution  of  Burgers  Equation 
Second  Order  Crank  Nicolson 


V 

N 

DX 

2 

\>DT/DX 

MAX 

ERROR 

L2 

ERROR 

L2  rate 

50 

.20 

6.25 

.633 

* 

io-3 

.125 

* 

io“2 

.500 

100 

.10 

25.00 

.158 

k 

10-3 

.311 

* 

io-3 

2.007 

200 

.05 

100.00 

.395 

k 

io~4 

.778 

* 

io-4 

1.999 

50 

.20 

3.125 

.303 

k 

io-2 

.442 

* 

io-2 

2.020 

.250 

100 

.10 

12.50 

.747 

k 

10  3 

.109 

* 

10  2 

200 

.05 

50.00 

.186 

k 

io-3 

.273 

* 

io-3 

1.997 

50 

.20 

1.5625 

.128 

k 

io-1 

.131 

* 

io-1 

.125 

100 

.10 

6.25 

.303 

k 

io-2 

.314 

k 

io-2 

2.061 

200 

.05 

25.00 

.749 

k 

io-3 

.775 

k 

io-3 

2.019 

50 

.20 

.775 

.694 

k 

io-1 

.473 

k 

io-1 

.062 

100 

.10 

3.100 

.130 

k 

io-1 

.940 

k 

io-2 

2.331 

2.069 

200 

.05 

12.40 

.308 

k 

io-2 

.224 

k 

10  2 

100 

.10 

1.55 

.694 

k 

io-1 

.334 

k 

io-1 

.031 

200 

.05 

6.20 

.130 

k 

io-2 

.665 

k 

io-2 

2.328 
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Tabic  6.5 

Steady  State  Solution  of  Burgers  Equation 
OCI  Crank  Nicolson  & Lees 
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Table  6.6 

Steady  State  Solution  of  Burgers  Equation 
Comparison  of  U Profiles 

v = .500 


X 

EXACT 

U 

OCT  - C-N 

2ND  ORDER  C-N 
N = 200 

N = 10 

N = 20 

N = 100 

-5.00 

.993307 

.993307 

.993307 

.993307 

.993307 

c 

o 

<r 

1 

.982014 

.982042 

.982015 

.982014 

.982021 

-3.00 

.952574 

.952545 

.952589 

.952574 

.952595 

-2.00 

.880797 

.881716 

.880850 

.880797 

.880833 

-1.00 

.731059 

.732380 

.731138 

.731059 

.731094 

- .40 

.598688 

.598688 

.598705 

- .20 

.549834 

.549834 

.549842 

-0.00 

.500000 

.500000 

. 500000 

.500000 

.500000 
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Table  6.7 


Steady  State  Solution  of  burgers  Equation 
Comparison  of  U Prof  lies 

v = .031 


X 

EXACT 

U 

OCT 

CN 

2ND  ORDER  CN  , 

l 

N » 200  1 

N = 60 

N *=  100 

-1.200 

1.000000 

1.000000 

1 . oooooo  j 

-1.167 

1.000000 

. 999997 

-1.000 

1 . oooooo 

. 999990 

1.000000 

1. OOOOOO  j 

- .833 

. 099999 

.999963 

; - .son 

.999998 

.999995 

1.000000 

- .667 

. 999979 

.999894 

- .600 

. 999937 

.999903 

1 . oooooo 

- .500 

.999086 

.999651 

- .400 

.998425 

.998091 

.999981 

- .333 

.995397 

.998843 

- .200 

.961794 

.962779 

.994937 

- .167 

.936325 

.996115 

0.000 

. 500000 

.500000 

. 500000 

. 500000 

.167 

.063675 

.003885 

j .200 

.038206 

.037221 

.005062 

.333 

.004603 

.001157 

.400 

.001575 

.001909 

.000019 

. 500 

.000314 

.000349 

.600 

.000063 

.000097 

.000000 

.667 

.000021 

.000106 

.800 

.000002 

.000005 

.oooooo 

.833 

.000001 

.000032 

1 000 

.000000 

.000010 

.000000 

.000000 

1.167 

.000000 

.000003 

1.200 

.000000 

. oooooo 


.OOOOOO 

2 


NSWC/WOL/TP  77-29 


Table  6.8 

Two  Dimensional  Parabolic  Equation 
OCI  - Crank-Ni colson 

u = a(x,y,t)  u + b(x,y,t)  u + c(x,y,t)  u + d(x,y,t)  u 
t ■ xx  x yy  y 

u(x,y,t)  = exp  { (xd-1) (y+1) (t+1) } 

Domain  is  square  fi  = [1/2  1 x,y  1 1]  , Ax  = Ay  = h 


1 TIME 
STEPS 


5 

20 

80 


.0 
40 
L60 


MAX 

MAX 

RELATIVE 

RELATIVE 

h 

At 

L - ERROR 

L -RATE 

ERROR 

RATE 

-1 

.1 

3.235-03 

4.414 

1.544-04 

3.905 

.05 

.025 

1.517-04 

4.955 

1.031-05 

3.977 

.025 

.00625 

4.890-06 

6.549-07 

.1 

3.903-02 

4.364 

3.909-04 

3.933 

.025 

1.896-03 

4.909 

2.559-05 

3.982 

.00625 

6.311-05 

1.619-06 

3 


<1 


Blasius  Flow 
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Hamming  formula  used  at 
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Table  6.10 
Howarth  Flow 


■ 

METHOD 

N 

An 

c + 

sep 

RUNNING  TIME 
(SEC) 

Second  Order 

20 

.50 

.1252 

8.90 

Crank  Nicolson 

50 

.20 

.1209 

19.92 

100 

.10 

.1204 

38.41 

200 

.05 

.1202 

71.58 

2x2  Block 

20 

.50 

.1202* 

19.38 

Tridiagonal 
Crank  Nicolson 

50 

.20 

.1202 

47.47 

100 

.10 

.1202 

94.05 

OCI 

20 

.50 

.1188* 

13.74 

Crank  Nicolson 

50 

.20 

.1199* 

35.81 

100 

.10 

.1201* 

73.63 

OCI 

20 

.50 

.1186** 

11.66 

Lees 

50 

.20 

.1200 

27.55 

100 

.10 

.1200 

54.21 

AC  = 10“A  , e***  = 10' 


+ Point  where  the  shear  changed  sign 
* Calculation  didn't  converge  and  was  terminated 
**  Oscillations  occurred  and  calculation  was  terminated 
***  Iteration  parameter 

r 

| 


* 
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APPENDIX  A 

The  operator  compact  implicit  formulas  are  derived  here  for  uniform  and 
nonuniform  grids,  with  their  associated  truncation  errors. 

G 1 ven 

(A-l)  Lu  = au  + bu  , 

XX  X 

a linear  relationship  between  u and  the  spatial  operator,  Lu,  at  x.  is 
sought  in  the  form 

fA-21  r u_  + r^  u^  + r+  u+  = q L(u)  + q L(u)  , + q+  I.(u)  + 

where  as  shorthand  notation  the  subscripts  -,  0,  + are  used  for  j-1,  j and  j+1 
respectively,  and  the  j dependence  of  the  coefficients  is  not  indicated,  see  (2.12). 


The  function  values  u_  and  u+  and  the  spatial  operators  L(u)_  and  L(u)+  can  be 
obtained  through  Taylor ' s series  expansion  about  the  point  j. 

«->  h+  ^ ^3)  ♦ ^ + + £ ^6,+-  • • 
«-»>  ■?’ ♦ £ ^ ^ * £ -T  - ^5)  * ^ 


(A-3c)  L(u) 


(A- 3d) 


a_  i/2)  + b_  u^!)  = b_  „<U  + (a_  - h_  b_)u^2) 


- h (a-  “ T b-)u^)  + TT  (a-  ' Y1  b-)  uo 

3 r (a-  - r h-)u!>5)  + It  (a-  - f b-)u 


(6) 

0 


T , , (2)  , (1)  . (1)  . , . . , , (2) 

L(u)+  = a+  u+  + b+  u+  = b+  un  + (a+  + h+  b+)un 


(A) 

'0 


h+^+ + t "n” + it  (a+ + r b+)“ 


where  superscripts  in  parenthesis  indicate  derivatives  and  h+  = x^+j  ~ xj  ancI 
h * j.  Multiplying  equations  (A-3a)  - (A-3d)  by  n,  R,  y,  6,  respectively  and 


collecting  terms  the  following  relation  is  obtained 


A-l 
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(A-4)  otu+  + Pu_  + yL(u)_  + AL(u)+  = (a+B)u0  + B + A u^2)  + C u^3)  + D u^4) 

or 

(A-5)  au+  - (m+fi)u0  + Pu_  = -yL(u)_  - AL(u)+  + A u^2)  + R U(!j13  + C u^3)  + B 


+ Truncation  Error  , 

where  In  order  that  (A-5)  be  equal  to  (A-2) , 


B r nb+  - Rh_  + A b+  + Yb 


0 


2 2 
ah.  Ph 


A = 


(A-fi) 


3 3 

ah  Bh 

c = ~y, 3T-  + «h. 


D = 


ah  4 Ph  4 Ah  2 


(a+  + h+  b+)  + Y (a-  ' h-  b-)  = a0 
'+  (a+  + r h+) _ Yh-  (a-  ■ 2^  0 = 0 


Def Ine 

(A-7)  a = aD,  B = BB,  y = yB,  A = AB 
where  B Is  the  determinant  of  (A-6) 

B = |l2a_a+(h+3  + 4h+  h + 4h+  h_  + h ) 

- 2a+b_h_(3h+3  + 7h42  h_  + 5h+  h_2  + h_3) 

(A-8)  « 23 

+ 2a_b+h+(h+  + 5h+^  h_  + 7h+  h_  + 3h_  ) 

- h+h_b+h_(h+  + h_)3}  , 

A A * A 

then  the  variables  a,  R,  y and  A are  Riven  by 


(A-9) 


y = |l2h+  a+  a0  Ch+2  - h+  h_  - h_2) 

+ 2a+bQh+2h_(3h+  + 2h_)  + 2a0b+h+2(h+2  - h+h_  - 2h_2) 
+ h+3h_b0b+(h+  + h_)| 


A-2 
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(A-in)  S = J 1 2a^a_h_ (h_‘  - h+h_  - h+?)  + ^0b_h_? (?h+?  + h+h  - h_2) 
- 2a_b0h+h_2(2h+  + 3h_)  + b0b_h+h_\h4  + h_)| 


(A-l  1 ) 

Pb_ (h+  + 

h_)  = D 

2an  - h+y 

- fi 

2a+  + h+b+]  - y [2a_ 

- b_ 

(2h_  + h+)] 

(A-12) 

Ph+  (h+  + 

h_)  = n 

2a0  + h_hn_ 

- 6 

2a+  + b+(2h+  + h_)J 

2a_  - h_b_j 

Multiplying  thru  by  -0,  the  q's  and  r's  become 

(A-13)  q = y , q+  = 4 , q^  = - D 

- * 0 * * + 
r = -fi,  r = (a  + R) , r = -« 

such  that  the  operators  0 and  R are  given  in  the  form 


(A-l 4) 


o = 6s+-nT  + ys_ 

R = - aR+  + (a  + B)T  - fi  R 


where  R is  the  shift  operator. 


Using  the  relations  (A-7)  - (A-12)  the  truncation  error  given  by 


i oh. 


(A-l 5) 


Rh 


FT  d)  5! 


5! 


] ( oh+^  Bh  ^ 4h4^ 

D » ~6T  + TT  + ~ 4 


ir  (a+  + r b+)  - (a-  - r b-)}""5 

r + r b+) + tt  - r b-)| uo 


is  seen  to  be  third  order  accurate  for  small  h. 

For  a uniform  mesh,  h+  = h = h,  the  truncation  error  reduces  to 


(A-lfi) 


FT  “ 1800  a 


- 35  a_a+b^  + 4 af)a+b_ 
which  is  fourth  order  accurate. 


jr  {-  [9  aoa-a+]u 

Kl 


(6) 

0 


+ 4 a0a_b+ 


Ncte  that  in  equation  (2.12),  common  factors  in  the  q's  and  r's  have  been 
canceled  (involving  constant  h) , so  that  (2-12)  differs  from  (A-13)  by  a 

3 

multiplicative  constant,  2h  . 


A- 3 
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APPENDIX  R 


The  coefficients  of  (6.6)  namely 


(B-l) 


(U  )n  * 
x 0 


Vo  + "l°l  + H2U2  + 


C,  f + 

0 0 


Vl  + C2f2 


are  derived. 


Consider  the  compact  implicit  formulas 


(B-2a) 
(B-2b) 
( B- 2 c ) 


Fn  “V[2'5  <u2  - V 


S0  + 10S1 + s2  - rf  % - 2li  + V 

h 


s0  + 1S1  + ?2  ■ iT  <F2  - V 


and  the  differential  equation  at  points  j=0,l,2  expressed  as 
(B-3)  a.S  + b ,F , = f.  j=0,l,2  . 

i i j j 1 

Equations  (B-2)  - (B-3)  form  a system  of  6 equations  in  9 unknowns,  and 
thus  Fq  can  be  determined  as  a function  of  uq> » U2  ’ ^1  an<^  ^2' 

The  coefficients  in  (B-l)  are  listed  below. 


(B-4) 


R-l 


I 
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The  truncation  error  is  given  by 


TRUNC 


(^)(^  - 10  9]  - Hh| 


where 


F = s il_  „(6)  = h^  (5)  F _ hi  (6) 

ES  5 30  ’ EF  30  ’ ET  20  u 

Eq.  (B-5)  can  be  specialized  for  constant  coefficients 


(B-6)  E = - (5  - + 77 

TRUNC  £)  180  L\  3 h 

a 


)“<6)+(w)"i 


In  the  case  of  time  dependent  problems  modifications  to  (B-l)  are  necessary. 
Consider  the  one  dimensional  parabolic  equation 


( B—  7 ) 


u = L(u)  - aS  + bF  = f 


The  first  derivative  at  an  end  point  at  time  level  (n+1)  in  the  form  of 
(B-8)  F"+1  = H;+1  Uj+1  + Hf 1 U^1  + H2n+1  uf 1 + 


„n+l  n+1  „n+l  Fn+1  , ~n+l  ^n-t  1 
G0  f0  + G1  fl  + C’2 


is  sought. 


Again,  as  before  use  the  compact  Implicit  formulas  (B-2) , but  with  u,F  and 
S evaluated  at  time  level  (n+1).  The  differential  equation  (B-7) , however, 
discretized  temporally  by  a Crank  Nicolson  scheme  to  yield 


l/?+1  - U° 


(B-9) 


[■r1  sfJ  + *] 

[*; + b"  « 5 5 


fn+!  + fn 


NSWC/WOL/TR  77-29 


Thus  f appearing  in  ( R—  8 ) is  the  spatial  operator  evaluated  at  (n+1)  and  is 
given  bv 


(B-10) 


rn+l  _ _n+l  n+1  , n+1  = 2_  .,n+l  _ r n 2_  m 

j ‘j  D.i  At  j 1 j At  ]> 


hence,  substituting  (B-10)  into  (B-8) , the  desired  relationship  is  obtained 

c1  • [«r ♦ h °r]  "S+1  + [»r 1 + it  <']  "T1 

(B-ll)  + [n"+1  + |j  Cf1]  if1  * Cj+1  [f”  + U"] 

+ <r  [r; + a!  “3 + Bf 1 [f2 + if  °;] 


The  local  spatial  truncation  error  remains  unchanged. 


